Pseudotime analysis of T-cells in developing thymus

library(Seurat)
library(conos)

Loading required package: igraph

Attaching package: ‘igraph’

The following object is masked from ‘package:scater’:

    normalize

The following objects are masked from ‘package:DelayedArray’:

    path, simplify

The following objects are masked from ‘package:dplyr’:

    as_data_frame, groups, union

The following objects are masked from ‘package:purrr’:

    compose, simplify

The following object is masked from ‘package:tidyr’:

    crossing

The following object is masked from ‘package:tibble’:

    as_data_frame

The following object is masked from ‘package:GenomicRanges’:

    union

The following object is masked from ‘package:IRanges’:

    union

The following object is masked from ‘package:S4Vectors’:

    union

The following objects are masked from ‘package:BiocGenerics’:

    normalize, path, union

The following objects are masked from ‘package:stats’:

    decompose, spectrum

The following object is masked from ‘package:base’:

    union

library(ggpubr)
library(tidyverse)
library(SingleCellExperiment)
# library(monocle3)
source("~/multiOmic_benchmark/utils.R")
source("~/multiOmic_benchmark/integrateBenchmark.R")

Based on the results of my benchmark, I set out to align expression and accessibility profiles from the F74 developing thymus dataset to detect changes in accessibility along pseudotime trajectories. While the benchmark was based on the task of label propagation, I here use the two most faithful methods (Seurat CCA and Conos) to achieve a common embedding of ATAC-seq and RNA-seq cells.

Load datasets.

rna.sce <- readRDS("~/my_data/F74_RNA_seurat_processed.RDS")
atac.sce <- readRDS("~/my_data/F74_ATAC_snapAtac_processed_bgmat.RDS")

Filter genes with zero variance

rna.gene.var <- as.matrix(counts(rna.sce)) %>% rowVars()
atac.gene.var <- as.matrix(counts(atac.sce)) %>% rowVars()

rna.sce <- rna.sce[which(rna.gene.var > 0),]
atac.sce <- atac.sce[which(atac.gene.var > 0),]

rna.sce; atac.sce

class: SingleCellExperiment 
dim: 24510 8321 
metadata(0):
assays(3): counts cpm logcounts
rownames(24510): RP11-34P13.3 RP11-34P13.7 ... AC233755.1
  AC240274.1
rowData names(0):
colnames(8321): AAACCTGAGTTCGATC_1 AAACCTGCAAGTTGTC_1 ...
  TTTGTCAAGCTGAACG_2 TTTGTCAGTATTAGCC_2
colData names(1): annotation
reducedDimNames(0):
spikeNames(0):
class: SingleCellExperiment 
dim: 31122 5793 
metadata(0):
assays(3): counts cpm logcounts
rownames(31122): A1BG A1BG-AS1 ... ZYX ZZEF1
rowData names(0):
colnames(5793): AAACGAAAGTGAACCG-1 AAACGAACATCGGCCA-1 ...
  TTTGTGTTCGATCGCG-1 TTTGTGTTCTGAGTAC-1
colData names(28): orig.ident nCount_ATAC ... nFeature_ACTIVITY
  ident
reducedDimNames(2): LSI UMAP
spikeNames(0):

Integration of T cells clusters

I re-run the integration based on the T cell subset. To select cells from the scATAC dataset, I take the SnapATAC clusters that best correspond to T-cells, based on label transfer.

tcells.sce.atac <- atac.sce[,which(as.numeric(atac.sce$seurat_clusters) %in% c(1:9))]

tcells.rna.ix <- which(rna.sce$annotation %in% c("DN","DP (Q)", "DP (P)", "SP (1)", "SP (2)"))
tcells.sce.rna <- rna.sce[,tcells.rna.ix]

tcells.sce.list <- list(RNA=tcells.sce.rna, ATAC=tcells.sce.atac)

Next, I select genes on which to perform integration. I take the union of the most variable features in the RNA dataset and the most covered features in the ATAC dataset

Remove cell cycle genes, that might interfere with pseudotime ordering

# integrate_features_ref <- hvg.rna 
# integrate_features_ref <- setdiff(integrate_features_ref, cell_cycle_genes) 
liger.obj <- scaleNotCenter(liger.obj, remove.missing = F)

Error in intI(i, n = x@Dim[1], dn[[1]], give.dn = FALSE) : 
  invalid character indexing

Visualize T cells in RNA dataset

Visualize markers

Visualize T cells in ATAC dataset: is the trajectory visible in the binary matrix?

Run CCA and Conos

Run CCA

coembed <- merge(x = seurat.list$RNA, y = seurat.list$ATAC)

coembed <- ScaleData(coembed, features = integrate_features_union, do.scale = FALSE)

Centering data matrix

  |                                                                                                          
  |                                                                                                    |   0%
  |                                                                                                          
  |==============                                                                                      |  14%
  |                                                                                                          
  |=============================                                                                       |  29%
  |                                                                                                          
  |===========================================                                                         |  43%
  |                                                                                                          
  |=========================================================                                           |  57%
  |                                                                                                          
  |=======================================================================                             |  71%
  |                                                                                                          
  |======================================================================================              |  86%
  |                                                                                                          
  |====================================================================================================| 100%

coembed <- RunPCA(coembed, features = integrate_features_union, verbose = FALSE)
coembed <- RunUMAP(coembed, dims = 1:30)

10:53:47 UMAP embedding parameters a = 0.9922 b = 1.112
10:53:47 Read 12078 rows and found 30 numeric columns
10:53:47 Using Annoy for neighbor search, n_neighbors = 30
10:53:47 Building Annoy index with metric = cosine, n_trees = 50
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
10:53:50 Writing NN index file to temp file /tmp/RtmplOAmA4/file511c460e48ef
10:53:50 Searching Annoy index using 1 thread, search_k = 3000
10:53:55 Annoy recall = 100%
10:53:59 Commencing smooth kNN distance calibration using 1 thread
10:54:03 Initializing from normalized Laplacian + noise
10:54:03 Commencing optimization for 200 epochs, with 562314 positive edges
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
10:54:17 Optimization finished

coembed <- AddMetaData(coembed, metadata = ifelse(colnames(coembed) %in% colnames(seurat.list[[reference]]), reference, query), col.name = "tech")

DimPlot(coembed, group.by = c('tech', "annotation"))

Run Conos

FeaturePlot(coembed, features = t.cell.markers$known.markers, split.by = "tech", slot = "data", cols = viridis::viridis(n=100))

Transfer labels on ATAC dataset

FeaturePlot(coembed, features = "prediction.score.max", cells = which(coembed$tech=="ATAC")) + scale_color_viridis_c()

Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale.

Run Pseudotime analysis

Identify cell of origin based on IGLL1 and CD34

AddMetaData(coembed, ifelse(colnames(coembed)==cell.oo, TRUE, FALSE), col.name = "iroot_cell")

An object of class Seurat 
55632 features across 12078 samples within 2 assays 
Active assay: RNA (24510 features)
 1 other assay present: ATAC
 2 dimensional reductions calculated: pca, umap

merged.sce <- SingleCellExperiment(list(counts=coembed@assays$RNA@counts, logcounts=coembed@assays$RNA@data), colData=coembed@meta.data[, c("annotation", "tech", "iroot_cell")],
                     reducedDims = map(coembed@reductions, ~ .x@cell.embeddings))

saveRDS(object = merged.sce, "~/my_data/Tcells_CCA_integration_20191127.RDS")
saveRDS(object = integrate_features_union, "~/my_data/intFeatures_Tcells_CCA_integration_20191127.RDS")

dpt <- read.csv('~/my_data/Tcells_CCA_integration_20191127_scanpy_dpt.csv') %>%
  select(X, dpt_pseudotime)

coembed <- AddMetaData(coembed, column_to_rownames(dpt, 'X'))

Check expression of markers along pseudotime

Bin ATAC cells by pseudotime

Fraction of accessible bins at each pseudotime bin

Motif analysis

Call peaks

peaks.ls = mclapply(seq(clusters.sel), function(i){
  print(paste("cluster", clusters.sel[i]))
  peaks = runMACS(
      obj=snap.out[which(snap.out@metaData$barcode %in% colnames(tcells.sce.atac)[tcells.sce.atac$seurat_clusters==clusters.sel[i]]),], 
      output.prefix=paste0("Tcells_F74_", gsub("cluster", clusters.sel)[i]),
      path.to.snaptools="/opt/conda/bin/snaptools",
      path.to.macs="/opt/conda/bin/macs2",
      gsize="hs", # mm, hs, etc
      buffer.size=500, 
      num.cores=3,
      macs.options="--nomodel --shift 100 --ext 200 --qval 5e-2 -B --SPMR",
      tmp.folder=tempdir()
 )
peaks
}, mc.cores=5)

all scheduled cores encountered errors in user code

Using chromVAR implementation in snapATAC

# snap.out = makeBinary(x.sp, "pmat")
snap.out@mmat = runChromVAR(
  obj=snap.out[which(snap.out@metaData$barcode %in% dpt.bin.cells),],
  input.mat="bmat",
  genome=BSgenome.Hsapiens.UCSC.hg38,
  # min.count=10,
  species="Homo sapiens"
);

Epoch: checking depedent packages ... 
Epoch: checking input parameters ... 
Epoch: creating chromVAR object ... 
Epoch: computing GC bias ... 
Epoch: getting JASPAR motifs ... 


There are 1 result in use. The connection will be released when they are closed